library(ggplot2)         #for data visualisation 
library(marginaleffects) #for plotting marginal effects
library(mclogit)         #for conditional logit models
library(MetBrewer)       #for colour palettes
library(patchwork)       #for combining plots
library(showtext)        #for displaying fonts in data viz
library(stargazer)       #for model tables in latex
library(survival)        #alternative package for conditional logistic regressions
library(tidyverse)       #for data-wrangling

##original-data-----
BrSenSwitch <- read.csv("BrSenSwitch.csv", sep = ";") %>%
               mutate(fusao = ifelse(is.na(fusao), 0, fusao)) %>%
               mutate(fissao = ifelse(is.na(fissao), 0, fissao)) %>%
               mutate(partido_novo = ifelse(is.na(partido_novo), 0, partido_novo)) %>%
               mutate(home_inverted = 1-home) %>%
               filter(partido_novo != 1,
                      fissao != 1) %>%
               mutate(to_be_elected = case_when(legislatura == 48 ~ 1,
                                                legislatura == 49 ~ 2,
                                                legislatura == 50 ~ 1,
                                                legislatura == 51 ~ 2,
                                                legislatura == 52 ~ 1,
                                                legislatura == 53 ~ 2,
                                                legislatura == 54 ~ 1,
                                                legislatura == 55 ~ 2,
                                                legislatura == 56 ~ 1))


double_renovation <- BrSenSwitch %>% filter(to_be_elected == 2) 
single_renovation <- BrSenSwitch %>% filter(to_be_elected == 1)

##merging different datasets----
#parties' policy positions based on roll-call votes
buoy <- BrSenSwitch %>%
        select(cardapio, legislatura, policy1) %>%
        unique() %>%
        drop_na() %>%
        #harmonising labels to make it easier to compare across datasets
        mutate(cardapio = case_when(cardapio == "PMDB" ~ "MDB",
                                    cardapio == "PFL" ~ "DEM",
                                    cardapio == "PPS" ~ "CID",
                                    cardapio == "CIDADANIA" ~ "CID",
                                    cardapio == "PP_2" ~ "PPold",
                                    cardapio == "PPR" ~ "PP",
                                    cardapio == "PRB" ~ "REP",
                                    cardapio == "REPUBLICANOS" ~ "REP",
                                    cardapio == "PCdoB" ~ "PCDOB",
                                    cardapio == "PR" ~ "PL",
                                           T ~ cardapio)) %>%
        mutate(let = paste0(cardapio, "_", legislatura))

#parties' policy positions based on elite survey (BLS)--
#This dataset can be downloaded from: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/6KVTUV
load("bls9_estimates_partiespresidents_long.RData")

nav_aid <- long.table %>%
           mutate(legislatura = case_when(year == 1990 ~ 48,
                                          year == 1993 ~ 49,
                                          year == 1997 ~ 50,
                                          year == 2001 ~ 51,
                                          year == 2005 ~ 52,
                                          year == 2009 ~ 53,
                                          year == 2013 ~ 54,
                                          year == 2017 ~ 55,
                                          year == 2021 ~ 56)) %>%
           mutate(let = paste0(party.or.pres,  "_", legislatura)) %>%
           select(let, ideo)

uniao <- nav_aid %>% 
         filter(let == "DEM_56" | let == "PSL_56") %>%
         summarise(let = "UNIÃO_56",
                   ideo = sum(ideo) / 2)

nav_aid <- nav_aid %>% rbind(uniao)

wtower <- buoy %>%
          left_join(nav_aid, by = "let") %>% 
          arrange(legislatura, policy1) %>%
          drop_na()


mdb_55 <- wtower %>%
          filter(let == "MDB_55") %>%
          mutate(policy1 = mean(policy1)) %>%
          unique()

wtower <- wtower %>%
          filter(let != "MDB_55") %>%
          rbind(mdb_55)
    
fogsinal <- wtower %>%
            pivot_longer(cols = c(policy1, ideo), names_to = "source", values_to = "policy")  %>%
            mutate(source = ifelse(source == "policy1", "Roll-Call\nVotes (PECs)", "BLS")) %>%
            mutate(legislatura = case_when(legislatura == 48 ~ "48th Legislature",
                                           legislatura == 49 ~ "49th Legislature",
                                           legislatura == 50 ~ "50th Legislature",
                                           legislatura == 51 ~ "51st Legislature",
                                           legislatura == 52 ~ "52nd Legislature",
                                           legislatura == 53 ~ "53rd Legislature",
                                           legislatura == 55 ~ "55th Legislature",
                                           legislatura == 56 ~ "56th Legislature"))           

##fonts----
#These fonts can be obtained here: https://fonts.google.com/specimen/Barlow+Semi+Condensed
font_add(family = "regular", "Barlow Semi Condensed-Regular.ttf")
font_add(family = "bold", "BarlowSemiCondensed-SemiBold.ttf")
showtext_auto() 

##Dataviz----
fogsinal %>% ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(linewidth = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 7, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties")) +
             theme(
                ##Titles
                plot.title=element_text(family="bold", hjust=0.5, vjust = 0.5, size=25, color="black"),
                plot.title.position = "plot",
                plot.subtitle = element_text(family="bold", size=15, hjust=0.5, color="black"),
                plot.caption= element_text(family="bold", size=15, color="black", hjust=0.5),
                plot.caption.position = "plot",
                ##Background
                #panel.border=element_blank(),
                panel.grid.major.y = element_blank(),
                panel.grid.major.x = element_blank(),
                #panel.grid.minor = element_line(colour = "grey"),
                plot.background = element_rect(fill = "white", linetype = 'blank'),
                panel.background = element_rect(fill = "white"),
                #plot.margin=ggplot2::margin(0.5, 0.5, 0.5, 0.5, "in"),
                #Axes
                axis.ticks.length=unit(0.15, "cm"),
                axis.ticks = element_blank(),
                axis.line = element_blank(),
                axis.title = element_text(size=30, family="bold", color="black"),
                axis.text.y = element_text(size=15, family="regular", color="black"),
                axis.text.x = element_text(size=15, family="regular", color="black"),
                ##Legend
                legend.position = "right",
                legend.title = element_text(size=20, family="regular"),
                legend.text = element_text(size=15, family="regular"),
                legend.background = element_rect(fill = "white"),
                #Strip
                #strip.background = element_rect(fill = "#fdf6ee"),
                strip.text = element_text(size=20, family="bold", hjust=0.5, vjust=0.5)) 


##Deviant cases----
fogsinal %>% filter(legislatura == "48th Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

fogsinal %>% filter(legislatura == "49th Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

th49 <- fogsinal %>% 
        filter(legislatura == "49th Legislature") %>% 
        arrange(source, policy) %>%
        ##deviant cases
        filter(let %in% c("PSB_49", "PRN_49"))
               
fogsinal %>% filter(legislatura == "50th Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

fogsinal %>% filter(legislatura == "51st Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

th51 <- fogsinal %>% 
        filter(legislatura == "51st Legislature") %>% 
        arrange(source, policy) %>%
        ##deviant cases
        filter(let %in% c("CID_51"))

fogsinal %>% filter(legislatura == "52nd Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

th52 <- fogsinal %>% 
        filter(legislatura == "52nd Legislature") %>% 
        arrange(source, policy) %>%
        ##deviant cases
        filter(let %in% c("PDT_52", "MDB_52", "PL_52", "PTB_52"))


fogsinal %>% filter(legislatura == "53rd Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

th53 <- fogsinal %>% 
        filter(legislatura == "53rd Legislature") %>% 
        arrange(source, policy) %>%
        ##deviant cases
        filter(let %in% c("PSOL_53", "MDB_53", "PP_53", "PL_53", "PTB_53"))

fogsinal %>% filter(legislatura == "55th Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))

th55 <- fogsinal %>% 
        filter(legislatura == "55th Legislature") %>% 
        arrange(source, policy) %>%
        ##deviant cases
        filter(let %in% c("PDT_55", "CID_55"))

fogsinal %>% filter(legislatura == "56th Legislature") %>%
             ggplot(aes(x = source, y = policy, group = cardapio))+
             geom_line(size = 1.5)+
             geom_point(aes(fill = cardapio), shape = 21, size = 5.5, stroke = 1.1)+
             xlab("") + ylab("") +
             facet_wrap(~legislatura) +
             scale_fill_manual(values = met.brewer("Signac", 24)) + 
             scale_color_manual(values = met.brewer("Signac", 24)) +
             guides(fill = guide_legend(title = "Parties"))


th56 <- fogsinal %>% 
        filter(legislatura == "56th Legislature") %>% 
        arrange(source, policy) %>%
        ##deviant cases
        filter(let %in% c("UNIÃO_56"))

th <- th49 %>% 
      rbind(th51, th52, th53, th55, th56) %>%
      filter(source != "BLS") %>%
      mutate(let = ifelse(let == "CID_51", "PPS_51", let),
             let = ifelse(let == "MDB_52", "PMDB_52", let),
             let = ifelse(let == "PL_53", "PR_53", let),
             let = ifelse(let == "MDB_53", "PMDB_53", let),
             let = ifelse(let == "CID_55", "PPS_55", let),
             let = ifelse(let == "CID_55", "PPS_55", let)) %>%
      select(let) %>%
      unique() %>%
      pull()

##Discarding deviant cases from the original dataset----
`%notin%` <- Negate(`%in%`)

BrSenSwitch_ideo_check <- BrSenSwitch %>%
                          mutate(let = paste0(cardapio, "_", legislatura)) %>%
                          filter(let %notin% th)

BrSenSwitch_ideo_check_double_renovation <- BrSenSwitch_ideo_check %>% filter(to_be_elected == 2) 
BrSenSwitch_ideo_check_single_renovation <- BrSenSwitch_ideo_check %>% filter(to_be_elected == 1)

##Robustness check----
model_a <- mclogit(
  y|choice ~ office + dif1 + party_str + 
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check)

model_b <- mclogit(
  y|choice ~ alt_swps + dif1 + party_str + 
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check)

model_c <- mclogit(
  y|choice ~ alt_swps + dif1 + party_str*elec_comp + 
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check)

model_d <- mclogit(
  y|choice ~ alt_swps + dif1 + party_str + 
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check_single_renovation)

model_e <- mclogit(
  y|choice ~ alt_swps + dif1 + party_str + party_str*elec_comp +
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check_single_renovation)

model_f <- mclogit(
  y|choice ~ alt_swps + dif1 + party_str + 
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check_double_renovation)

model_g <- mclogit(
  y|choice ~ alt_swps + dif1 + party_str*elec_comp + 
    home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted,
  data = BrSenSwitch_ideo_check_double_renovation)

##models with clogit----
model1 = clogit(y ~ office + dif1 + party_str + 
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check)

model2 = clogit(y ~ alt_swps + dif1 + party_str + 
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check)

model3 = clogit(y ~ alt_swps + dif1 + party_str*elec_comp + 
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check)

model4 = clogit(y ~ alt_swps + dif1 + party_str + 
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check_single_renovation)


model5 = clogit(y ~ alt_swps + dif1 + party_str + party_str*elec_comp +
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check_single_renovation)


model6 = clogit(y ~ alt_swps + dif1 + party_str + 
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check_double_renovation)

model7 = clogit(y ~ alt_swps + dif1 + party_str*elec_comp + 
                  home + coligacao_governador + suplente*home_inverted + legislaturas_participadas*home_inverted + pop_rural*home_inverted +
                  strata(choice), data=BrSenSwitch_ideo_check_double_renovation)

##stargazer----
stargazer(model1, model2, model3, model4, model5, model6, model7, type = "text")

stargazer(model_a, model_b, model_c, model_d, model_e, model_f, model_g, 
          column.labels   = c("Full Dataset", "M = 1", "M = 2"),
          column.separate = c(3, 2, 2),
          no.space=T, order = c(1,2,3,4,7,5,6,8,9, 10), digits = 2, 
          covariate.labels = c("Presidential Coalition", "Weighted Participation in Coalition", 
                               "Ideological Distance", "Party Strength", "Party Strength * Electoral Competition",
                               "Home Party", "Governor's Coalition", "Surrogate", "Number of Legislatures",
                               "Rural Population"),
          ##input from clogit models for more detailed infmation
          add.lines=list(c("N","211,106","178,859","178,859", "95,913", "95,913", "82,946", "82,946"),
                         c("R$^{2}", "0.330", "0.336", "0.336", "0.324", "0.324", "0.349", "0.349"),
                         c("Max. Possible R2", "0.335", "0.340", "0.340",  "0.329", "0.329", "0.352", "0.352"),
                         c("Log Likelihood", "-754.970", "-536.005",  "-535.732", "-346.029", "-346.024", "-179.968", "-179.866")))